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ABSTRACT 

We present a new method for calculating linear cosmic microwave background 

(CMB) anisotropy spectra based on integration over sources along the photon 

^ 1 past light cone. In this approach the temperature anisotropy is written as a 

m ' 
O 

geometrical term is given by radial eigenf unctions which do not depend on the 
particular cosmological model. The source term can be expressed in terms of 
photon, baryon and metric perturbations, all of which can be calculated using a 
small number of differential equations. This split clearly separates between the 
dynamical and geometrical effects on the CMB anisotropies. More importantly, 



Of 



2 ' it allows to significantly reduce the computational time compared to standard 

methods. This is achieved because the source term, which depends on the model 
and is generally the most time consuming part of calculation, is a slowly varying 
function of wavelength and needs to be evaluated only in a small number of 
points. The geometrical term, which oscillates much more rapidly than the 
source term, does not depend on the particular model and can be precomputed 
in advance. Standard methods that do not separate the two terms and require 
a much higher number of evaluations. The new method leads to about two 
orders of magnitude reduction in CPU time when compared to standard 
methods and typically requires a few minutes on a workstation for a single 
model. The method should be especially useful for accurate determinations of 
cosmological parameters from CMB anisotropy and polarization measurements 
that will become possible with the next generation of experiments. A programm 
implementing this method can be obtained from the authors. 



Subject headings: cosmology: cosmic microwave background, cosmology: 
large-scale structure of the universe, gravitation, cosmology: dark matter 



-2- 



1. Introduction 



The field of cosmic microwave background (CMB) anisotropies has seen a rapid 
development since its first detection by the COBE satellite only a few years ago. There 
are now several reported experimental results that are detecting anisotropies on degree 
angular scales (see [Scott, Silk fc White 1995| and [Bond 19961 f° r a recent review), which 



together with a few upper limits on smaller angular scales already give interesting limits 
on cosmological models. With the development of the new generation of experiments 
now being proposed one hopes to accurately map the CMB sky from arcminute scales to 
several degree scales. The amount of data thus provided would allow for an unprecedented 
accuracy in the determination of cosmological parameters. As theoretical modelling shows 
( |Bond et al. T99g [Hu, Sugiyama fc Silk 19951 , [Jungman et al. 1995| , [Seljak 1994p , CMB 



anisotropies are sensitive to most of the cosmological parameters and have a distinctive 
advantage over other cosmological observations in that they probe the universe in the linear 
regime. This avoids the complications caused by physical processes in the nonlinear regime 
and allows to use powerful statistical techniques to search over the parameter space for the 
best cosmological model (see e.g. Jungman et al. 1995|) . A large stumbling block in this 



program at present is the speed of theoretical model calculations, which are still too slow 
to allow for a rapid search over the parameter space. This limitation was partially removed 
by the development of appoximation methods ([Hu fc Sugiyama 19954 b; |Seljak 1994|) , 



which can give fast predictions of CMB anisotropy with a 10% accuracy. However, these 
approximations are not accurate enough to exploit the complete amount of information 
that will be present in the future CMB observations. This is especially true for some of the 
more extreme cosmological models, where simple approximations break down and lead to 
systematic inaccuracies in the results. Obviously, it would be useful to have a fast method 
that would not be based on any approximations and would lead to accurate results for 
any cosmological model. The purpose of this paper is to present a new method of CMB 
calculation that satisfies these requirements. 

Theoretical calculations of the CMB anisotropies are based on linear theory of 
cosmological perturbations, developed first by Lifshitz (1946) and applied to the CMB 
anisotropies by Peebles and Yu (1970). In this early calculation only photons and baryons 
were included, but later workers extended the calculations to include dark matter (Bond & 
Efstathiou 1984, 1987; Vittorio & Silk 1984), curvature (Wilson & Silk 1981; Sugiyama & 
Gouda 1992; White & Scott 1995), gravity waves or tensor modes (Crittenden et al. 1993) 
and massive neutrinos (Bond & Szalay 1983; Ma & Bertschinger 1995; Dodelson, Gates & 
Stebbins 1995). Most of these and more recent calculations (e.g. Holtzmann 1989; Stompor 
1994; Sugiyama 1995) solve for each Fourier mode of temperature anisotropy Ay(fc) by 
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expanding it in Legendre series up to some desired l max and then numerically evolve this 
system of equations in time from the radiation dominated epoch until today. Typically this 
means evolving a system of several thousand coupled differential equations in time, a slow 
process even for the present day computers. In addition, because each multipole moment is 
a rapidly oscillating function one has to densely sample it in values of k with typical number 
of evaluations of the order of l ma x- Even the fastest codes at present require several hours 
of CPU time for each theoretical model (Sugiyama 1995), while some numerically more 
accurate ones (e.g. Pode fc Bertschinger 1995|) require more like tens or hundreds of hours. 

In this paper we explore a different approach to compute CMB anisotropies based on 
integration of the sources over the photon past light cone. The method is a generalization 
of approximate method developed by one of the authors (Seljak 1994). It differs from it 
in that it is exact, in the sense that it can achieve arbitrary precision within the limits of 
linear perturbation theory. By rewriting the system of equations in the integral form one 
can separate between the geometrical and dynamical contributions to the anisotropies. The 
former do not depend on the model and need to be computed only once, while the latter 
contain all the information on the particular model and can be computed with a small 
system of equations. Solving for CMB anisotropies using this integral form greatly reduces 
the required computational time. The outline of the paper is as follows: in §2 we present 
the basic system of equations that needs to be solved both in the standard and in the 
integral method. In §3 we present in some detail a practical implementation of the integral 
method, highlighting the computational differences between it and the standard Boltzmann 
method. We conclude in §4 by discussing possible applications where the new method can 
be particularly useful. 



2. Method 

In this section we first present the standard system of equations that needs to be solved 
for temperature anisotropies, which is based on solving Boltzmann equation using Legendre 
expansion of photon distribution function. This part follows closely the existing literature 
(e.g. Ma & Bertschinger 1995, Bond 1996 and references therein) and only the final results 
are given. We do not discuss the technical details of the standard Boltzmann method, 
except where our approach differs significantly from it. In the second part of the section we 
present the integral solution of photon distribution, which is the basis of our method. In 
this paper we restrict the analysis to a spatially flat universe. 
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2.1. Boltzmann, Einstein and Fluid Equations 



The temperature anisotropy at position x in the direction n is denoted with Aj(f,n). 
In principle it depends both on the direction and on the frequency, but because spectral 
distortions are only introduced at the second order the frequency dependence can in the 
lowest order be integrated out. Anisotropy At{x, n) can be expanded in terms of Fourier 
modes Ay(A;,n), which in linear perturbation theory evolve independently of one another. 
Assuming perturbations are axially-symmetric around k, we may further Legendre expand 
the anisotropy in the angle fx — k ■ ft/k, 

A T (k,n)=J2( 2l + 1 )H) lA TiPi(n), (1) 
i 

where Pi(ji) is the Legendre polynomial of order I and A?i is the associated multipole 
moment. A similar decomposition also applies to the amplitude of polarization anisotropy 
Ap(k, n) ( |Bond fc Efstathiou 1984] , Critenden et al. 1993, Kosowsky 1995, Zaldarriaga & 



Harari 1995). 

Evolution of the temperature anisotropy is governed by the Boltzmann equation 
( Peebles fc Yu 1970| , |Wilson fc Silk 1981j , Pond fc Efstathiou 1984Q . Its collisionless part is 



given by the time component of the geodesic equation, which depends on the metric. Here 
we will use the metric in the longitudinal gauge ( |Bardeen 1980| ; Bertschinger 1996), which 



is similar to the gauge-invariant formalism (Kodama & Sasaki 1984) and gives expressions 
that are most similar to their Newtonian counterparts. The choice of gauge is purely a 
matter of convenience and in some cases (e.g isocurvature models) other gauge choices such 
as synchronous gauge are computationally advantageous over the longitudinal gauge (see 
e.g. Pode fc Bertschinger 1995] , Pond 1996]) . In the longitudinal gauge the perturbations 
are specified with two scalar potentials <fi and ip and a gauge-invariant tensor perturbation 
h (we will ignore vector perturbations in this paper, as they most likely have a negligible 
contribution to CMB anisotropy). The corresponding temperature and polarization 

( S) ( S) (T) (T) 

anisotropies are denoted as A^ , A P for scalar and Aj. , Ap for tensor components. In 
linear perturbation theory the scalar and tensor perturbations evolve independently and 
the total power is given by the sum of the two contributions. 

The collisional part of the photon Boltzmann equation is determined by the Thomson 
scattering term. After angular and momentum integration the Boltzmann evolution 
equations for scalar perturbations can be written as (Pond fc Efstathiou 1987]), 



A^ s) + ikfiA^ = <t>- ikfiip + k{-A { T } + A^ + ifiv b + ^P 2 (/u)II} 
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Here the derivatives are taken with respect to the conformal time r and v b is the velocity 
of baryons. Differential optical depth for Thomson scattering is denoted as k — an e x e <7T, 
where a(r) is the expansion factor normalized to unity today, n e is the electron density, x e 
is the ionization fraction and ot is the Thomson cross section. The total optical depth at 
time t is obtained by integrating k, k{t) = k(r)dr. A useful variable is the visibility 
function g{r) = kexp(—n). Its peak defines the epoch of recombination, when the dominant 
contribution to the CMB anisotropies arises. 

Expanding the temperature anisotropy in multipole moments one finds the following 
hierarchy of coupled differential equations ( |Wilson fc Silk 1981| , [Bond fc Efstathiou 1984| , 
|Ma fc Bertschinger 1995|) , 
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where <5y is the Kronecker symbol. A similar system of equations without the Thomson 
scattering terms and polarization also applies for massless neutrinos. For the massive 
neutrinos the system of equations is more complicated, because the momentum dependence 
cannot be integrated out of the expressions (see e.g. |Ma fc Bertschinger 1995| ). 



Baryons and cold dark matter can be approximated as fluids and their evolution can be 
obtained from the local conservation of energy-momentum tensor. This gives the equations 
for cold dark matter density 5 C and its velocity v c , 



-kv c + 



■ - v c + kip . 
a 



(4) 



For baryons one has additional terms in the Euler's equation caused by Thomson scattering 
and pressure, 



Vb 



-kv b + 
a 
a 



v b + c 2 s k5 b + ^k(3A$ - v b ) + kip . 



(5) 
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Here c s is the baryon sound speed and p 7 , pt, are the mean photon and baryon densities, 
respectively 

Finally, the evolution of scalar metric perturbations is given by Einstein equations, 
which couple the sources and the metric perturbations. Only two equations are needed to 
specify the evolution. Here we choose them to be the energy and momentum constraint 
equations, 

k 2 4> + 3- U + -rj>) = -AnGa 2 5p 
a \ a J 

:2 i : , a „/. \ i .2 



k (0 + -t/jj = 4irGa 2 5f, (6) 

where 5p and Sf are the total density and momentum density perturbations, respectively. 
They are obtained by summing over the contributions from all species, 5p = J2i$Pi, 
5f = J2i $fii $Pi — Pi$i an d dfi — (Pi + Pi)Vi, where pi and pi are the mean density and 
pressure of the i-ih species. 

For tensor perturbations the Boltzmann equation is given by (Crittenden et al. 1993), 

A^ T) + ikpA ( p = -h - k(A ( p - *) , 
AJP + ikpA ( p = -k(A ( p + *) , 

(7) 



* = 



1a (t) + — A (T) + — A (T) - -A (T) + — A (T) - — A (T) 



10 ' u 35 " z 210 ^ 5 ™ 35 ^ z 210 



The only external source is that of the tensor metric perturbation which evolves according 
to the Einstein equations as 

h + 2-h + k 2 h = 0. (8) 

a 

We ignored the source term on the right-hand side of equation above (caused by neutrino 
and photon anisotropic stress), as it is always negligible compared to the terms on the 
left-hand side. 

— * 

To obtain the temperature anisotropy for a given mode k one has to start at early 
time in the radiation dominated epoch with initial conditions of the appropriate type 
(e.g. isentropic or isocurvature) and evolve the system of equations until the present. The 
anisotropy spectrum is then obtained by integrating over the initial power spectrum of the 
metric perturbation P^(k), 

C\ S) = (4vr) 2 j k 2 dkP4k)\A i T S l \k,r = r )| 2 . (9) 

Analogous expression holds for the polarization spectrum and for the tensor spectrum 
(where the initial power spectrum P^(k) has to be replaced by the initial tensor power 
spectrum Ph{k)). 
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The spectrum Ci is related to the angular correlation function, 

1 oo 

C(9) = (A{n 1 )A(n 2 )) nvH2=cose = -r-E( 2/ + l)C,P,(«>s0). (10) 

To test a model on a given angular scale 9 one has to solve for A77 up to I ~ 1/0. If one 
is interested in small angular scales this leads to a large system of differential equations 
to be evolved in time and the computational time becomes long. For a typical spectrum 
with I max ~ 1000 one has to evolve a system of 3000 differential equations (for photon 
and neutrino anisotropy and photon polarization) until the present epoch. Moreover, the 
solutions are rapidly oscillating functions of time, so the integration has to proceed in small 
time increments to achieve the required accuracy on the final values. 



2.2. Integral solution 

Instead of solving the coupled system of differential equations (|3|) one may formally 
integrate equations § along the photon past light cone to obtain (e.g. Zaldarriaga & Harari 
1995), 

a£? } = r dre ik ^ T - To) e- K {ke- K [A TO + i/iv b + -P 2 (^M + - ik^} 
Jo 2 

A { p ] = -- r dre ik ^ T - To) e- K k[l - P 2 (/i)]n. (11) 
2 Jo 

Expressions above can be further modified by eliminating the angle \x in the integrand 
through the integration by parts. The boundary terms can be dropped, because they vanish 
as t — ► and are unobservable for r = r (i.e. only the monopole term is affected). This 
way each time a given term is multiplied by a fi, it is replaced by its time derivative. This 
manipulation leads to the following expression, 

A^ P = r dre ikfl{T - To) S { T S p(k,T) 
Jo 
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Some of the terms in the source function (r) are easily recognizable. The first two 
contributions in the first term are the intrinsic anisotropy and gravitational potential 
contributions from the last-scattering surface, while the third contribution is part of the 
velocity term, the other part being the k~ 1 gvb term in the second row. These terms make 
a dominant contribution to the anisotropy in the standard recombination models. The 
first term in the second row, e~ K (0 + if>), is the so-called integrated Sachs- Wolfe term and 
is important after recombination. It is especially important if matter-radiation equality 
occurs close to the recombination or in fi ma ttcr 7^ 1 models. In both cases gravitational 
potential decays with time, which leads to an enhancement of anisotropies on large 
angular scales. Finally we have the terms caused by photon polarization and anisotropic 
Thomson scattering, which contribute to II. These terms affect the anisotropy spectra at 
the 10% level and are important for accurate model predictions. Moreover, they are the 
sources for photon polarization. Equation ([12]) is a generalization of the tight-coupling 
and instantaneous recombination approximation (Seljak 1994) and reduces to it in the 
limit where the visibility function is a delta-function and II can be neglected. In that 
approximation one only needs to evaluate the sources at recombination and then free 
stream them to obtain the anisotropy today. In the more general case presented here one 
has to perform an additional integration over time, which includes the contributions arising 
during and after the recombination. Moreover, because the tight-coupling approximation 
is breaking down at the time of recombination, both polarization and photon anisotropic 
stress are being generated and II makes a non-negligible contribution to the anisotropy. 
For exact calculations one has to use ([12]), which properly includes all the terms that are 
relevant in the linear perturbation theory. 

To solve for the angular power spectrum one has to expand the plane wave e lk ^ T ~ T °^ 
in terms of the radial and angular eigenf unctions (spherical Bessel functions and Legendre 
polynomials, respectively), perform the ensemble average [] and integrate over the angular 
variable /x. This leads (|9|), where the multipole moment at present time A^ P y(k,T = r ) is 
given by the following expression, 



where ji(x) is the spherical Bessel function. Note that while angular eigenfunctions 



^n performing ensemble average we assume that only the amplitude and not the phase of a given mode 
evolves in time. While this is valid in linear theory for most models of structure formation, it may not be 
correct in some versions of models with topological defects (Albrecht et al. 1995). 




(k, t = r ) = / S^Uk, T)ji[k(r - r)]dr, 
Jo 



(13) 
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( s) 

Sj> p, which does not depend on the multipole moment I and a geometrical term ji, which 
does not depend on the particular cosmological model. The latter thus only needs to be 
computed once and can be stored for subsequent calculations. The source term is the same 
for all multipole moments and only depends on a small number of contributors in ( O) 
(gravitational potentials, baryon velocity and photon moments up to I = 4). By specifying 
the source term as a function of time one can compute the corresponding spectrum of 
anisotropies. Equation flilf ) is formally an integral system of equations, because / < 4 
moments appear on both sides of equations. To solve for these moments it is best to use 
the equations in their differential form (0), instead of the integral form above. Once the 
moments that enter into the source function are computed one can solve for the higher 
moments by performing the integration in ([13]) (see section 3 for more details). 



The solution for the tensor modes can similarly be written as an integral over the 
source term and the tensor spherical eigenfunctions \\. The latter are related to the 
spherical Bessel functions ( [Abbott fc Schaeffer 1986]) , 



This gives 



(l + 2)\ 3l (kr) 



A (T) 



\ 2(1-2)1 (kr) 2 ' 

rdrSPp(k,T) X l k (T -T), (15) 
Jo 



where from equation [7] follows 

4 T > = -he~ K + g^> S [ P = -gV. (16) 



Equations (p^2| - |l6D are the main equations of this paper and form the basis of the line of sight 



integration method of computing CMB anisotropies. In the next section we will discuss in 
more detail the computational advantages of this formulation of Boltzmann equation and 
its implementation. 



3. Calculational Techniques 



In the previous section we presented the expressions needed for the implementation of 
the line of sight integration method. As shown in ([13]) and (|i~5[) one needs to integrate over 
time the source term at time r multiplied with the sperical Bessel function evaluated at 
k(r — r). The latter does not depend on the model and can be precomputed in advance. 
Fast algorithms exist which can compute spherical Bessel functions on a grid in k and I in 
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short amount of time (e.g. Press et al. 1993). The grid is then stored on a disk and used 
for all the subsequent calculations. This leaves us with the task of accurately calculating 
the source term, which determines the CMB spectrum for a given model. Below we discuss 
some of the calculational techniques needed for the implementation of the method. We 
especially highlight the differences between this approach and the standard Boltzmann 
integration approach. Our goal is to develop a method which is accurate to 1% in C\ up 
to I ~ 1000 over the whole range of cosmo logical parameters of interest. These include 
models with varying amount of dark matter, baryonic matter, Hubble constant, vacuum 
energy, neutrino mass, shape of initial spectrum of perturbations, reionization and tensor 
modes. The choice of accuracy is based on estimates of observational accuracies that will 
be achiavable in the next generation of experiments and also on the theoretical limitations 
of model predictions (e.g. cosmic variance, second order effects etc.). Most of the figures 
where we discuss the choice of parameters are calculated for the standard CDM model. 
This model is a reasonable choice in the sense that it is a model which exhibits most of 
the physical effects in realistic models, including acoustic oscillations, early-time integrated 
Sachs- Wolfe effect and Silk damping. One has to be careful however not to tune the 
parameters based on a single model. We compared our results with results from other 
groups (Bode & Bertschinger 1995; Sugiyama 1995) for a number of different models. We 
find a better than 1% agreement with these calculations over most of the parameter space 
of models. The computational parameters we recommend below are based on this more 
detailed comparison and are typically more stringent than what one would find based on 
the comparison with the standard CDM model only. 



3.1. Number of coupled differential equations 

In the standard Boltzmann method the photon distribution function is expanded to a 
high / max (^) and typically one has to solve a coupled system of several thousand differential 
equations. In the integral method one evaluates the source terms S(k, r) as a function of 
time (0), flnp and one only requires the knowledge of photon multipole moments up to 
I = 4, plus the metric perturbations and baryon velocity. This greatly reduces the number 
of coupled differential equations that are needed to be solved. For an accurate evaluation 
of the lowest multipoles in the integral method one has to extend the hierarchy somewhat 
beyond / = 4, because the lower multipole moments are coupled to the higher multipoles 
(0). Because power is only being transferred from lower to higher I it suffices to keep a 
few moments to achieve a high numerical accuracy of I < 5 moments. One has to be 
careful however to avoid unwanted reflections of the power being transferred from low I 
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Fig. 1. — CMB spectra produced by varying the number of evolved photon multipole 
moments, together with the relative error (in %) compared to the exact case. While using 
L = 5 produces up to 2% error, using L = 7 gives results almost identical to the exact case. 



to high I, which occur for example if a simple cut-off in the hierarchy is imposed. This 
can be achieved by modifying the boundary condition for the last term in the hierarchy 
using the free streaming approximation ( [Ma fc Bertschinger 1995| , [Hu et al. 1995|) . In the 
absence of scattering (the so-called free streaming regime), the recurrence relation among 
the photon multipoles in equation (|3|) becomes the generator of spherical Bessel functions. 
One can therefore use a different recurrence relation among the spherical Bessel functions 
to approximate the last term in the hierarchy without reference to the higher terms. The 
same approximation can also be used for polarization and neutrino hierarchies. This type 
of closure scheme works extremely well and only a few multipoles beyond / = 4 are needed 
for an accurate calculation of the source term. This is shown in figure [I], where a relative 
error in the spectrum is plotted for several choices of maximal number of photon multipoles. 
We choose to end the photon hierarchy (both anisotropy and polarization) at / 7 = 8 and 
massless neutrino at l v = 7, which results in an error lower than 0.1% compared to the 
exact case. Instead of a few thousand coupled differential equations we therefore evolve 
about 35 equations and the integration time is correspondingly reduced. 



3.2. Sampling of CMB multipoles 
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1 

Fig. 2. — Relative error between the exact and interpolated spectrum, where every 20th, 
50th or 70th multipole is calculated. The maximal error for the three approximations is less 
than 0.2, 0.4 and 1.2%, respectively. The rms deviation from the exact spectrum is further 
improved by finer sampling, because the interpolated spectra are exact in the sampled points. 
For the sampling in every 50th multipole the rms error is 0.1%. 



In the standard Boltzmann integration method one solves for the whole photon 
hierarchy (^) and the resultant A; is automaticaly obtained for each / up to some / max . The 
CMB spectra are however very smooth (see figure except for the lowest /, where the 
discrete nature of the spectrum becomes important. This means that the spectrum need not 
be sampled for each I and instead it suffices to sparsely sample the spectrum in a number of 
points and interpolate between them. Figure |2] shows the result of such interpolation with 
cubic splines (see e.g. |Press et al. 1992|) when every 20th, 50th or 70th I is sampled beyond 
/ = 100 with an increasingly denser sampling towards small I, so that each I is sampled 
below / = 10. While sampling of every 70th I results in maximal error of 1%, sampling in 
every 20th or 50th / gives errors below 0.2 and 0.4%, respectively. We choose to compute 
every 50th C\ beyond / = 100 in addition to 15 / modes below / = 100, so that a total of 
45 / modes are calculated up to / max = 1500. This gives a typical (rms) error of 0.1%, with 
excursions of up to 0.4%. The number of integrals in equation (Jffy is thus reduced by 10-50 
and the computational time needed for the integrals becomes comparable or smaller than 
the time needed to solve for the system of differential equations. 
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3.3. Free streaming 



After recombination and in the absence of a time changing gravitational potential the 
source function often becomes negligible. This is the so called free streaming regime, where 
the photons are freely propagating through the universe. Most of the standard Boltzmann 
codes use a special free streaming algorithm to map the anisotropies from a given epoch Tf s 
into anisotropies today (|Bond fc Efstathiou 1984]) . In the line of sight integration method 
the free streaming regime is only a special case where S(k, r) = after some time Tf s . Thus 
one can stop the integration at the time Tf S beyond which the sources are not important 
and there is no need for a separate algorithm to evolve the anisotropies until today. For 
example, if one assumes that only the first term in ([12]) is important one only needs to 
integrate over the source where the visibility function g appreciably differs from 0. In 
the absence of reionization this restricts the time integration to a narrow interval during 
recombination around z ~ 1100. Although most of the contributions to the anisotropies 
come from this epoch, time dependent gravitational potential (and, to a smaller extent, 
other source terms in equation (0)) make a nonnegligible contribution to the anisotropies 
even after recombination. As mentioned earlier, this is especially important in f2 mat t er 7^ 1 
models and in models with low Q matt evh 2 . In the first case the gravitational potential is 
decaying at late times, while in the second class of models the matter-radiation equality 
during which gravitational potential changes in time is pushed to a lower redshift. Even 
in standard CDM model (f2 matter = l,h = 0.5) gravitational potential is still significantly 
changing in time at moderately low redshifts of z ~ 100 ( |Hu et al. 1995| ). Similarly one 
cannot use free streaming in the models with late reionization, where the visibility function 
is nonvanishing at low redshifts. We choose to integrate until the present time for most of 
the models, except for the models with fi matte r = 1? where we stop the integration at z = 10. 
In this case the computational time is reduced significantly (typically 50%) compared to 
that of evolving the equations until the present time. 



3.4. Integration over time 



For each Fourier mode k the source term is integrated over time r QTB"!). The sampling 
in time need not be uniform, because the dominant contribution arises from the epoch of 
recombination around z ~ 1100, the width of which is determined by the visibility function 
g and is rather narrow in look-back time for standard recombination scenarios. During 
this epoch the sources acoustically oscillate on a time scale of c s fc _1 , so that the longest 
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wavelength modes are the slowest to vary. For short wavelengths the rate of sampling 
should therefore be higher. Even for long wavelengths the source function should still be 
sampled in several points across the last-scattering surface. This is because the terms 
in (|T^) depend on the derivatives of the visibility function. If the visibility function g is 
narrow then its derivative will also be narrow and will sharply change sign at the peak of g. 
Its integration will lead to numerical roundoff errors if not properly sampled, even though 
positive and negative contributions nearly cancel out when integrated over time and make 
only a small contribution to the integral. Figure § shows the error in integration caused by 
sampling this epoch with 10, 20 or 40 points. Based on comparison with several models 
we choose to sample the recombination epoch with 40 points, which results in very small 
(~ 0.1%) errors. After this epoch the main contribution to the anisotropies arises from the 
integrated Sachs- Wolfe term. This is typically a slowly changing function and it is sufficient 
to sample the entire range in time until the present in 40 points. The exceptions here are 
models with reionization, where the visibility function becomes non-negligible again and a 
new last-scattering surface is created. In this case a more accurate sampling of the source 
is also needed at lower redshifts. 



3.5. Integration over wavenumbers 

The main computational cost of standard CMB calculations is solving the coupled 
system of differential equations. The number of fc-modes for which the system is solved 
is the main factor that determines the speed of the method. For results accurate to / max 
one has to sample the wavenumbers upto a maximun value fc max — ^max/ r o- 111 the line of 
sight integration method solving the coupled system of differential equations still dominates 
the computational time (although for each mode the time is significantly shorter than in 
the standard Boltzmann method because of a smaller system of equations). It is therefore 
instructive to compare the number of k evaluations needed in each of the methods to 
achieve a given accuracy in the final spectrum. 

In the standard Boltzmann method one solves for A^pJ(k) directly, so this quantity 
must be sampled densely enough for accurate integration. Figure [|a shows A^}(k) for 
/ = 150 in a standard CDM model. One can see that it is a rapidly oscillating function 
with a frequency k ~ Tq 1 . Each oscillation needs to be sampled in at least a few points to 
assure an accurate integration. To obtain a smooth CMB spectrum one typically requires 
6 points over one period, implying 2Z max k-mode evaluations. This number can be reduced 
somewhat by filtering out the sampling noise in the spectrum ( |Hu et al. 1995| ), but even in 



this case one requires at least 1-2 points per each period or / max /2 /c-mode evaluations. 
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Fig. 3. — Error in the spectrum caused by insufficient temporal sampling of the source term. 
Inaccurate sampling of the source during recombination leads to numerical errors, which can 
reach the level of 1% if the source is sampled in only 10 points across the recombination 
epoch. Finer sampling in time gives much smaller errors for this model. Comparisons with 
other models indicate that sampling in 40 points is needed for accurate integration. 



To understand the nature of these rapid oscillations in A T l(k) we will consider 
wavelenghts larger than the width of the last scattering surface. In this case the Bessel 
function in (|13|) can be pulled out of the integral as ji(kr ) because the time at which 
recombination occurs, when the dominant contribution to A^(k) is created, is much 
smaller than r and kAr 1 (At is the interval of time for which the visibility function 
differs appreciably from zero). So the final A^i(k) is approximately the product of 
ji(kro) and integrated over time, if the finite width of the last scattering surface and 
contributions after recombination can be ignored. 

Figure f|b shows the source term integrated over time and the Bessel function 
ji(kr ). It shows that the high frequency oscillations in A^ z (/c) seen in figure |]a are caused 
by the oscillation of the spherical Bessel functions, while the oscillations of the source 
term have a much longer period in k. The different periods of the two oscillations can be 
understood using the tight coupling approximation (|Hu et al. 1995| , [Seljak 19941 ). Prior and 
during recombination photons are coupled to the baryons and the two oscillate together 
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Fig. 4.— In (a) A^ ] l50 (k) is plotted function of wavevector k. In (b) A T ^ 50 (k) i 

( s) 

decomposed into the source term integrated over time and the spherical Bessel function 
ji5o(kr ). The high frequency oscillations of A^\ 50 (k) are caused by oscillations of the 
spherical Bessel function ji5o(^To), whereas the source term varies much more slowly. This 
allows one to reduce the number of k evaluations in the line of sight integration method, 
because only the source term needs to be sampled. 



(S) 



is 



with a typical acoustic timescale t s ~ r rec /\^3 ~ tq/\/3z icc ~ r /50. The frequency of 
acoustic oscillations k ~ is therefore 50 times higher than the frequency of oscillations 
in spherical Bessel functions, which oscillate as Tq 

Because an accurate sampling of the source term requires only a few points over 
each acoustic oscillation, the total number of k evaluations in the integral method can be 
significantly reduced compared to the standard methods. Typically a few dozen evaluations 
are needed over the entire range of k, compared to about 500 evaluations in the standard 
method when a noise filtering technique is used and 2000 otherwise (for / max ~ 1000). 
Once the source term is evaluated at these points one can interpolate it at points with 
preevaluated spherical Bessel functions, which can be much more densely sampled at no 
additional computational cost. The end result is the same accuracy as in the standard 
method, provided that the source is sampled in sufficient number of points. Figure |5| shows 
the relative error in the CMB spectrum for the cases where the source term is calculated 
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Fig. 5. — Error in the spectrum caused by insufficient /c-mode sampling of the source term. 
Sampling the source with 40 points up to k = 2/ max leads to 1% errors, while with 60 or 80 
points the maximal error decreases to 0.2%. Comparisons with other models indicate that 
sampling in 60 points is sufficient for accurate integration. 

in 40, 60 and 80 points between and kr = 3000 (for Z max = 1500). While using 40 
points results in up to 1% errors, using 60 points decreases the maximum error to below 
0.2% for this model. In general it suffices to use / max /30 k modes, which is at least an 
order of magnitude smaller than in the standard methods. Note that with this method 
there is no need to filter the spectrum to reduce the sampling noise, because the latter is 
mainly caused by insufficient sampling of the spherical Bessel functions, which are easy 
to precompute. The additional operations needed for a higher sampling (summation and 
source interpolation) do not significantly affect the overall computational time. Moreover, 
if each C\ is accurately calculated they can be sparsely sampled and interpolated (section 
3.2), this would not be possible if they had a significant noise component added to them. 



4. Conclusions 



In this paper we presented a new method for accurate calculations of CMB anisotropy 
and polarization spectra. The method is not based on any approximations and is an 
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alternative to the standard Boltzmann calculations, which are based on solving large 
numbers of differential equations. The approach proposed here uses a hybrid integro- 
differential approach in solving the same system of equations. By rewriting the Boltzmann 
equations in the integral form the solution for the photon anisotropy spectrum can be 
written as an integral over a source and a geometrical term. The first is determined by a 
small number of contributors to the photon equations of motion and the second is given by 
the radial eigenfunctions, which do not depend on the particular cosmological model, but 
only on the geometry of space. 

One advantage of the split between geometrical and dynamical terms is that it 
clarifies their different contributions to the final spectrum. A good example of this is 
the temperature anisotropy in the non-flat universe, which can be be written using a 
similar decomposition, except that spherical Bessel functions have to be replaced with 
their appropriate generalization QAbbott fc Schacffcr 1986 ). This will be discussed in more 



detail in a future publication, here we simply remark that replacing radial eigenfunctions 
in a non-flat space with their flat space counterpart (keeping comoving angular distance 
to the LSS unchanged) is only approximate and does not become exact even in the large I 
(small angle) limit. The geometry of the universe leaves its signature in the CMB spectra 
in a rather nontrivial way and does not lead only to a simple rescaling of the spectrum by 
fi ma ( tor (|Jungman et al. 1995|). 



The main advantage of our line of sight integration method is its speed and accuracy. 
For a given set of parameters it is two orders of magnitude faster than the standard 
Boltzmann methods, while preserving the same accuracy. We compared our results with 
the results by Sugiyama (1995) and by Bode & Bertschinger (1995) and in both cases the 
agreement was better than 1% up to a very high I for all of the models we compared to. 

The method is useful for fast and accurate normalizations of density power spectra 
from CMB measurements, which for a given model require the CMB anisotropy spectrum 
and matter transfer function, both of which are provided by the output of the method. 
Speed and accuracy are even more important for accurate determination of cosmological 
parameters from CMB measurements. In such applications one wants to perform a 
search over a large parameter space, which typically requires calculating the spectra of a 
several thousand models (e.g. |Jungman et al. 1995| ). One feasible way to do so is to use 



approximation methods mentioned in the introduction. These can be made extremely fast, 
but at a cost of sacrificing the accuracy. While several percent accuracy is sufficient for 
analyzing the present day experiments, it will not satisfy the requirements for the future 
all-sky surveys of microwave sky. Provided that foreground contributions can be succesfully 
filtered out (see Tegmark & Efstathiou 1995 for a recent discussion) one can hope for 
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accuracies on the spectrum close to the cosmic variance limit, which for a broad band 
averages can indeed reach below 1% at Z > 100. It is at this stage that fast and accurate 
CMB calculations such as the one presented in this paper will become crucial and might 
enable one to determine many cosmological parameters with an unprecedented accuracy. 

We would like to thank Ed Bertschinger for encouraging this work and providing 
helpful comments. This work was partially supported by grant NASA NAG5-2816. 
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